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Abstract 

A finite-size technique is employed to compute the normalization constant 
.Za of the isovector axial current in lattice QCD. The calculation is carried 
out in the quenched approximation for values of the bare gauge coupling g 
ranging from to 1. In the lattice action and the lattice expression for the 
axial current we include the counterterms required for 0(a) improvement, 
with non-perturbatively determined coefficients. With little additional work 
the normalization constant Zy of the improved isospin current is also obtained. 



February 2008 



1. Introduction 



In lattice QCD with Wilson quarks the conservation of the isovector axial 
current is violated by lattice effects. As a consequence a finite renormalization 
of the current is required to ensure that the chiral Ward identities assume 
their canonical form [1,2]. It is evidently important to compute the associated 
normalization constant Za, since it contributes directly to physical matrix 
elements such as the pion decay constant F^. 

In perturbation theory Za has been worked out to one-loop order for vari- 
ous lattice actions and lattice definitions of the axial current [3-7] . The results 
may be used to calculate Za at the couplings of interest, but since these are 
not small in general it is difficult to say how reliable the numbers are that one 
obtains. A non-perturbative determination of the normalization constant is 
clearly preferable. Two different strategies to perform such a calculation have 
been pursued. In the first case Za is fixed by requiring certain chiral Ward 
identities between correlation functions of the axial and vector currents to be 
satisfied on the lattice [2,8-10]. The correlation functions are then evaluated 
through numerical simulation. The other proposition is to compute matrix 
elements of the axial current between quark states and to determine the nor- 
malization of the current by matching the numerical results with renormalized 
perturbation theory at large momentum transfers [11]. 

Our principal aim in the present paper is to calculate Za in the on-shell 
0(a) improved lattice theory. The significance of improvement in this context 
has previously been stressed in refs. [12-14]. Here we employ the improved 
action and the improved axial current with non-perturbatively determined 
0(a) counterterms [15-18]. All calculations are carried out in the quenched 
approximation. We use the Ward identity method and combine it with a finite- 
size technique based on the Schrodinger functional. This allows us to set the 
quark mass to zero (or to values very close to zero) and to determine Za at 
all bare couplings go between and 1. Contact with perturbation theory can 
thus be made. 

For the definition of the Schrodinger functional and the 0(a) improved 
theory the reader is referred to ref. [16]. The notations introduced there are 
taken over completely without further notice. In sect. 2 we briefly recall the 
euclidean Ward identities associated with the chiral symmetry of QCD in the 
continuum limit. We then define the isospin vector and axial vector currents 
in the on-shell 0(a) improved lattice theory and derive the normalization con- 
ditions that will be used to compute the associated normalization constants 
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(sect. 3). As discussed in sect. 4 a careful interpretation of the role played 
by the current normalization conditions is required on the lattice, because the 
chiral Ward identities are only valid up to cutoff effects of order a 2 . The cal- 
culation of the isospin vector and axial vector current normalization constants 
through numerical simulations is described in sects. 5 and 6. The paper ends 
with a few concluding remarks and a technical appendix, where an essentially 
rigorous proof of the crucial Ward identity is given. 



2. Euclidean current algebra 

We first consider the theory in the continuum limit and proceed formally, 
i.e. without paying attention to the proper definition of the correlation func- 
tions that occur. The boundary conditions on the quark and gluon fields do 
not matter in this section. We assume that there is an isospin doublet of 
quarks with mass m and study the associated chiral symmetry of the theory. 

The isospin vector and axial vector variations of the quark and anti-quark 
field are defined by 

5W(x) = |r<V(x), 5^{x) = -#r)ir a , (2.1) 

8M(x) = lT Q 7 5^(x), 5 a 3(x) = #r) 7 5±r a , (2.2) 

where T a denotes a Pauli matrix acting on the flavour indices of the quark 
field. The definition extends to arbitrary expressions O by treating 5y and 
6% as first order differential operators. In particular, for the variations of the 
isospin vector and axial vector currents, 

V£(x) = #r) 7At ±T<V(x), Al(x) = #r) 7At75 ±T a V0r), (2.3) 

one obtains 

SvV>(x) = -ie abc VC(x), 6 a A V»(x) = ^6 abc ^(*), (2.4) 

<5«4(x) = -ie abc A;(x), SlA^x) = -ie abc V^(x). (2.5) 
The currents thus form a closed algebra under these variations. 
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The Ward identities associated with the chiral symmetry of the action 
are derived by performing local infinitesimal symmetry transformations of the 
quark and anti-quark fields in the euclidean functional integral. We choose 
to write the identities in an integrated form which is quite intuitive and will 
prove useful later on when we discuss the lattice theory 

Let R be a space-time region with smooth boundary dR. Suppose O int 
and O ex t are polynomials in the basic fields localized in the interior and exterior 
of R respectively. The general vector current Ward identity then reads 

/ da^x) (V«(x)O int O ext ) = - ((<5»0 int ) Oext) , (2.6) 

JdR 

while for the axial current one obtains 

/ da^(x) (Al(x)O int O ext ) = - ((5 a A O int ) O ext ) 

JdR 

+ 2m [ d 4 x {P a (x)O int O ext ) . (2.7) 

J R 

The integration measure da^(x) points along the outward normal to the sur- 
face dR and the pseudo-scalar density P a {x) is defined by 

P a (x)=^) 75 ir>(x). (2.8) 

The left-hand sides of eqs. (2.6) and (2.7) may be interpreted as matrix el- 
ements of the charge operators associated with the currents. This is made 
particularly clear if we choose R to be the space-time volume between two 
equal-time hyper-planes. 

For illustration we set the quark mass to zero and choose O mt to be one of 
the currents at some point y in the interior of R. An example of the resulting 
identities then is 

/ da,(x) (A*(x)Al(y)O ext ) = ie abc (V„ c (y)O ext ) , (2.9) 

JdR 

and three more such relations, corresponding to eqs. (2.4) and (2.5), may 
be obtained. In this way the current algebra, which one usually sets up in 
Minkowski space in an operator language, is recovered in the euclidean domain. 
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3. Currents in lattice QCD 



As already mentioned in sect. 1 we choose Wilson's formulation of lattice 
QCD and include all 0(a) correction terms that are required for on-shell im- 
provement. The lattice action etc. is exactly as in ref. [16]. We assume that 
the coefficients multiplying the 0(a) counterterms have been adjusted so that 
the residual cutoff effects are of order a 2 . 

3.1 Improved currents 

On the lattice the bare currents, and A^, are defined through the local 
expressions (2.3). The unrenormalized on-shell 0(a) improved currents are 
then given by [15,16] 

(V.Y^VZ + cyalidl + d^ (3.1) 

(AO" =Al + c A a\(dt + d,)P a . (3.2) 

For the anti-symmetric tensor field T£ v we can take 

T;(x)=#(x)a F ir>(x), (3.3) 

while the axial density P a is again defined through eq. (2.8). 

The renormalization of lattice QCD is particularly transparent if a mass- 
independent renormalization scheme is employed. As discussed in sect. 3 of 
ref. [16] the renormalized improved currents in such schemes are given by 

(V K )l = Z v (l + &vam q )(V02, (3.4) 
(A R )« = Z A (1 + 6 A am q )(A0^. (3.5) 

The renormalization constants Zy and Z A are functions of the modified bare 
coupling 

9o =5 , o( 1 + & g° m q) ) ( 3 - 6 ) 

while by, 6a and 6 g depend on go and should be adjusted so as to cancel 
any remaining cutoff effects of order am q . We shall not need to know these 
coefficients here, because the quark mass will be set to zero for the calculation 
of Zy and Z A . 
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3.2 Normalization condition for the vector current 

Although the isospin symmetry of the continuum theory is preserved on the 
lattice, the improved vector current introduced above is only conserved up to 
cutoff effects of order a 2 . Its normalization is hence not naturally given and 
we must impose a normalization condition to fix Zy. Our aim in the following 
is to derive such a condition by studying the action of the renormalized isospin 
charge on states with definite isospin quantum numbers. 

The matrix elements that we shall consider are constructed in the frame- 
work of the Schrodinger functional (see ref. [16] for details). We use the 
boundary field products 

O a = a 6 ^C(u)7 5 ir a C(v), (3.7) 

u.v 

0' a = a 6 ^C'(u)75|T«C'(v), (3.8) 

u,v 

to create initial and final states that transform according to the vector repre- 
sentation of the exact isospin symmetry. The correlation function 

x 

can then be interpreted as a matrix element of the renormalized isospin charge 
between such states. The charge generates an infinitesimal isospin rotation (if 
properly normalized) and after some algebra one finds that the correlation 
function must be equal to 

h = -^(O la O a ). (3.10) 

Strictly speaking this argumentation is only correct in the continuum theory. 
We may however conclude that 

/vVo) = /i + 0(a 2 ), (3.11) 

since the lattice correlation functions approach the continuum limit with a rate 
proportional to a 2 . Note that we do not need to include the renormalization 
factors for the boundary quark fields here because they cancel in eq. (3.11). 
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The 0(a) counterterm appearing in the definition (3.1) of the improved 
vector current does not contribute to the correlation function /y (xo)- So if 
we introduce the analogous correlation function for the bare current, 



Mx ) = — ^ abc {O' a V b (x)O c ), 



(3.12) 



X 



it follows from eq. (3.11) that 



Z v (l + b v am q )f v (x ) = /i + 0(a 2 ). 



(3.13) 



By evaluating the correlation functions /i and fv(xo) through numerical sim- 
ulation one is thus able to compute the normalization factor Zy(l + byam^). 
In particular, to calculate Zy it suffices to consider the theory at vanishing 
quark mass. 

3. 3 Normalization condition for the axial current 

To derive a normalization condition for Za, we set the quark mass to zero 
from the beginning. Our starting point is the Ward identity (2.9) which we 
now write in the form 



One may be hesitant to make use of this relation, since it has been deduced 
in a formal manner. A general argument, presented in appendix A, however 
shows that such worries are not justified. There is little doubt that eq. (3.14) 
is a true property of the theory in the continuum limit and so may be used to 
fix the normalization of the axial current on the lattice. 

We now pass to the lattice theory and assume Schrodinger functional 
boundary conditions as before. A convenient choice of the region R is the 
space-time volume between the hyper-planes at xq = yo ± t. From eq. (3.14) 
and 0(a) improvement we then expect that 




da^x) e abc (A«(x)Al(y)O ext ) = 2i (V„ c (y)O ext ) . 



(3.14) 



a: 



3 Y,^ (U R ) a (y + t,x) - (A R r (y - i,x 



)}(A R ) b (y)O ext 



> 



X 




(3.15) 
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Fig. 1. Quark diagrams contributing to fi (diagram a) and /v(xo) 
(diagrams b and c). Filled (open) circles represent the creation (annihila- 
tion) of a quark at the boundaries of the lattice. The squares indicate the 
vector current insertions. 

It has been important here that the fields in the correlation functions are 
localized at non-zero distances from each other. Since the theory is only on- 
shell improved, one would otherwise not be able to say that the error term is 
of order a 2 (cf. sect. 2 of ref. [16]). 

After summing over the spatial components of y, and using the fact that 
the axial charge is conserved at zero quark mass (up to corrections of order 
a 2 ), eq. (3.15) becomes 

a^e abc ((A R T (x)(A R ) b (y)O ext ) = 
x,y 

a 3 ^i((y R )§(y)a x t) + 0(a 2 ), (3.16) 
y 

where xq = yo + t. We now choose the field product O ext so that the function 
/v (z/o) introduced previously appears on the right-hand side of eq. (3.16). The 
normalization condition for the vector current then allows us to replace the 
correlation function by f\. Explicitly, we define 

/X A (x ,y ) = -^^e^e^^'^^OSW^O^y)^) (3-17) 

x,y 

and conclude from the above that 

Zlfl A (x ,y )=f 1 + O(a 2 ) (3.18) 

for all times xo > yo between and T. The normalization constant Za can 
thus be determined by computing the correlation functions f\ and f AA (x ,y ) 
at zero quark mass. 
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Fig. 2. Quark diagrams contributing to /| A (xo, yo)- The squares 
indicate the axial current insertions at time xq (upper square) and yo 
(lower square). 

3.4 Disconnected diagrams and the strange quark 

As usual the integration over the quark and anti-quark fields in the functional 
integral is carried out analytically. One is then left with an integration over 
all gauge fields, the correlation functions fi etc. being given by a set of quark 
diagrams that correspond to the possible Wick contractions of the quark and 
anti-quark fields. For Schrodinger functional boundary conditions the required 
two-point contractions have been worked out in detail in sect. 2 of ref. [17]. 

In the case of /1 and fy(xo) the assignment of the isospin quantum num- 
bers is such that no disconnected quark diagrams appear (see fig. 1). The cor- 
relation function /aa( x 0j 2/o) involves two current insertions and many more 
Wick contractions exist. Among the diagrams listed in fig. 2 there are also 
two disconnected ones (diagrams g and h) . The isospin factors associated with 
the diagrams e and f vanish so that they can be dropped immediately. 

We would now like to show that the disconnected diagrams do not con- 
tribute either. To this end we introduce a third quark, referred to as the 
strange quark, and replace the field products O la and O a in the definitions 
of f\ and /aa^Oj^o) by products of a strange and an isospin doublet quark 
boundary field. With appropriately contracted indices, the argumentation 



leading to eq. (3.18) then goes through unchanged. Note that the currents are 
the same as before and the calculated value of must hence come out to 
be the same up to terms of order a 2 . Since the strange quark does not cou- 
ple to the axial current, the number of possible Wick contractions is strongly 
reduced and only the diagrams a and c survive, the line without current in- 
sertions representing the strange quark propagator. If we interchange quarks 
and anti-quarks one obtains the diagrams b and d instead. The isospin fac- 
tors associated with the diagrams can be worked out straightforwardly and a 
comparison of the situation with and without strange quark then shows that 
the total contribution of the disconnected diagrams must vanish up to terms 
of order a 2 . 

In the course of the numerical computations described later in this paper, 
we have been able to verify that the disconnected diagrams indeed add up to 
zero within small statistical errors. We have thus decided to drop them and 
to extract from the connected part of f^(xo,yo). 



4. Lattice effects and current normalization 

The relations (3.13) and (3.18) determine the current normalization con- 
stants Zy and Za only up to cutoff effects of order a 2 . Depending on the 
choice of the lattice size, the boundary values of the gauge field and the other 
kinematical parameters that one has, slightly different results for Zy and Z& 
are hence obtained. One may try to assign a systematic error to the normal- 
ization constants by studying these variations in detail, but since there is no 
general rule as to which choices of the kinematical parameters are considered 
to be reasonable, such error estimates are bound to be rather subjective. 

In our opinion the better way to deal with the problem is to define the 
normalization constants through a particular normalization condition. The 
physical matrix elements of the renormalized currents that one is interested 
in must then be calculated for a range of lattice spacings so as to be able 
to extrapolate the data to the continuum limit. The results obtained in this 
way are guaranteed to be independent of the chosen normalization condition, 
because any differences in the normalization constants of order a 2 extrapolate 
to zero together with the cutoff effects associated with the matrix elements 
themselves. 
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In the following we shall adopt this point of view and the precise choices 
that we shall make are then not too important. Some care must be paid to 
ensure that the cutoff effects in matrix elements of the renormalized currents 
between low-energy states are not artificially enhanced through an inappropri- 
ate choice of the kinematical parameters in eqs. (3.13) and (3.18). In particu- 
lar, the external length scales (the time difference xq — yo and the lattice size, 
for example) should be sufficiently large compared to the lattice spacing, at all 
bare couplings considered. Perturbation theory can serve as a guide here and 
further confidence can be gained by studying the magnitude of the residual 
cutoff effects in various matrix elements of the renormalized currents through 
numerical simulations. 



5. Numerical evaluation of /i, fy and f AA 

In this section we present the details of the numerical calculations which 
we have performed to determine Za and Zy. We start in subsect. 5.1 by 
listing the expressions for /i, fy and f AA in terms of quark propagators. In 
subsect. 5.2 we briefly discuss the lattice action and algorithms used in the 
numerical simulation. We have closely followed the procedures outlined in 
ref. [18], which can be consulted for further information and any unexplained 
notations. 

5.1 Explicit form of f\, fy and f AA 

As already mentioned in subsect. 3.4, the mathematical expressions corre- 
sponding to the quark diagrams shown in figs. 1 and 2 are obtained by apply- 
ing Wick's theorem to the appropriate product of quark fields. There is only 
one diagram contributing to the correlation function f\ and one finds that 

f 1 = ±(tT{tfK}) G , (5.1) 

where the trace is over the Dirac and colour indices. The matrix K represents 
the quark propagation from the boundary at time to the boundary at time T. 
In terms of the solution H(x) of the lattice Dirac equation introduced in sect. 2 
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of ref. [18] it is given by 



The numerical calculation of H(x) is discussed in subsect. 3.2 of ref. [18]. 

In the case of the correlation function fy(xo), two diagrams contribute 
and one obtains 

/v(*o) = ^E< Re tr{^ t 7 5^'(x) t 757o^(x)}> G , (5.3) 



where H'(x) denotes the quark propagator from the boundary at time T to 
the point x in the interior of the space-time volume. H'(x) is defined through 
the Dirac equation 

(D + 5D + m )H'(x) = 0, < x < T, (5.4) 

and the boundary conditions 

P + H'(x)\ xo=0 = 0, P-H'(x)\ X0=T = P- (5.5) 

The numerical solution of these equations proceeds as in the analogous case 
of the matrix H(x) (cf. subsect. 3.2 of ref. [18]). 

Taking the definition (3.2) of the improved axial current into account, the 
correlation function / AA (x ,yo) m ay be expanded according to 

/aa = /aa + c A a [BSfpA + 4"/ap] + c 2 A a 2 dSd y f PP (5.6) 

with the obvious definitions of /aa, /ap ; /pa and /pp. The superscripts on 
the symmetric difference operator 

d, = \{dt + d,) (5.7) 

imply a differentiation with respect to x or y respectively. As discussed in 
subsect. 3.4 only the diagrams a-d in fig. 2 need to be evaluated for the com- 
putation of Za- The corresponding expressions for /aa. Zap, /pa and /pp 
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look very similar and we only give the result for /aa and diagram a, viz. 



{/AA(£o,yo)}diagram a = 

^3 £<tr {K^ 5 H'(x)^ 5l0l5 S(x,yh ^H(y)}) G . (5.8) 



The bulk quark propagator S(x,y) appearing in this formula is the inverse of 
the Dirac operator D + 5D + mo in the space of quark fields with vanishing 
boundary values. When evaluating /aa and the other correlation functions, 
the propagator itself is however not required. Instead one first calculates H(x) 
at all x and then solves the Dirac equation 

(D + dD + m )!.a 3 J2 S (^yhol 5 H(y)\=a- 1 S Xoyol ol 5 H(x) (5.9) 
^ y ' 

by applying the usual iterative methods. 
5.2 Details of the simulation 

As for the coefficients multiplying the 0(a) counterterms in the improved 
action and the improved axial current, we follow ref. [18] and set 

_ 1 - 0.656 g 2 - 0.152 gg - 0.054ff 6 
Csw " 1 - 0.922 gl ' (5 ' 1Uj 



rA ' uu07: ' 6g ° x La977^ (5 - n) 



for all couplings in the range < go < 1. These formulae have been ob- 
tained non-perturbatively by imposing some carefully chosen improvement 
conditions. We expect that an almost perfect cancellation of 0(a) effects in 
on-shell quantities is thus achieved. 

From now on we shall often quote values of (3 = Q/g^ and k = (8+2am ) _1 
instead of the bare coupling and mass. All our production runs have been per- 
formed on APE/Quadrics computers with 256 nodes. We used the same hybrid 
over-relaxation algorithm as described in subsect. 3.1 of ref. [18] to generate a 
representative ensemble of gauge field configurations. Subsequent evaluations 
of the correlation functions /i, fy and /aa were separated by 25 iterations 
of the algorithm. We have checked explicitly for the statistical independence 
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of our sample by dividing the full ensemble into bins, each containing a num- 
ber of individual "measurements" . The statistical errors were then monitored 
as the number of measurements per bin was increased. We did not observe 
any significant change of the errors for increasing bin size, which we take as 
evidence for the statistical independence of our sample. 

To invert the Dirac operator we employed the stabilized biconjugate gra- 
dient algorithm (BiCGstab) with even-odd preconditioning [19,20]. By com- 
paring our results with those from a set of Fortran-90 programs, we have 
verified not only the correct evaluation of correlation functions, but also that 
the rounding errors associated with the 32 bit arithmetic on the APE com- 
puter were completely negligible in our calculation. All statistical errors were 
estimated using the jackknife method. 



6. Computation of Zy and 

We now proceed to describe the non-perturbative determination of the 
normalization factors Zy and Za in the range < go < 1. Except for the 
tests mentioned in subsect. 6.4, the boundary values C and C of the gauge 
field and the angles 9k are set to zero throughout this section. 

6.1 Complete specification of the normalization conditions 

As discussed in sect. 4 we need to make a definite choice for the parameters 
on which the normalization conditions (3.13) and (3.18) depend. The quark 
mass is set to zero, as previously indicated, and the remaining parameters are 
then the times at which the currents are inserted and the lattice extensions T 
and L. We scale these parameters proportionally to L and eventually decided 
to take 

Z v /v(T/2) = A, T = 2L, (6.1) 
Z 2 A fl A (2T/3, T/3) = /i, T = 9L/4, (6.2) 

as the definite form of the normalization conditions. We still need to say, 
however, what precisely it means to set the quark mass to zero and how L/a 
is to be scaled with g . 
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The critical hopping parameter k c (i.e. the zero mass point) depends on 
the details of the lattice definition of the quark mass [16,18]. The ambiguity 
is just one of the sources of the order a 2 corrections in the normalization 
conditions for the currents and so is to be treated following the lines of sect. 4, 
i.e. we adopt any particular definition of the quark mass and use it to calculate 

K c . 

The definition that we have chosen is the same as the one previously em- 
ployed in sect. 7 of ref. [18]. The starting point is the unrenormalized current 
quark mass m{xo) which one extracts from the PCAC relation (eq. (5.2) of 
ref. [18]). We then set T = 2L and define k c to be the value of the hopping 
parameter k where 

2a 

m av = l rn(T/2 + t) (6.3) 

t=-2a 

vanishes. The average over the time coordinate xq is taken to reduce the 
statistical error on the calculated mass values [18]. Note that the critical 
hopping parameter so defined is slightly dependent on L/a. It is implicitly 
understood that at any given value of (3 one chooses a lattice size L/a and 
first computes k c and then evaluates the normalization conditions (6.1) and 
(6.2) at this value of k c and the same lattice size L/a. 

At this point the normalization constants Zy and are well-defined 
functions of go and L/a. The dependence on the lattice size is of order (a/L) 2 
in the continuum limit and we have verified that at the couplings of interest 
the change in the calculated values of the normalization constants is indeed 
small when increasing L/a from say 8 to 16. 

According to the discussion in sect. 4 we now need to make a definite 
choice of L/a. To ensure that the on-shell matrix elements of the renormalized 
improved currents approach the continuum limit with a rate proportional to 
a 2 , we must require that L remains fixed in physical units. Explicitly, we 
define L/a at all couplings go < 1 through 

L/a = 8 at g Q = 1, (6.4) 
L/vo = constant, (6.5) 

where r denotes a hadronic scale extracted from the force between heavy 
quarks [21]. Using recent lattice data for ro [22] one finds that L/a\/s=Q.2 — H, 
L/a\/3 = e.4 — 14 and L/a > 16 for (3 > 6.8. For practical reasons we did not 
perform the simulations at exactly these lattice sizes. In fact this is not really 
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Fig. 3. Simulation results for fi/ fv(xo) from a lattice of size 32 x 16 3 
at j3 = 6.4 and zero quark mass. Note that we are using a fine scale in this 
plot. The statistical fluctuations of the data points are on the level of a 
small fraction of a percent. 

required because the associated systematic errors can be estimated reliably 
and turn out to be small (details are given below). 

6. 2 Results for Zy and 

The statistical fluctuations of the numerators and denominators in the ratios 
fi/fv(xo) and /i/Jaa^o, Ho) are strongly correlated. The jackknife error 
estimation accounts for these correlations and the ratios are obtained with 
impressive statistical accuracy even if only a small ensemble of independent 
gauge field configurations is available. Fig. 3 shows the ratio fi/fv(xo) as a 
function of xq/ci from a typical run. One observes a clear signal with small 
statistical uncertainty and nearly no time-dependence within errors. The data 
points at different times xo are statistically decorrelated to such an extent 
that the signal-to-noise ratio can be enhanced by averaging the data in the 
range T/2 — 2a < xq < T/2 + 2a. Strictly speaking this should be taken as 
part of the definition of the normalization condition for Zy (in the sense of 
sect. 4), but we did not want to obscure the discussion in subsect. 6.1 with too 
many details and thus mention this item only now. In the case of the ratio 
/i//aa( x 0) Uo) n o averaging has been performed. 
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Table 1. The values of L/a and k c used in the normalization conditions 



L/a 




L/a 




6.0 8 


0.135046(16) 






6.2 8 


0.135692(6) 


12 


0.135705(12) 


6.4 8 


0.135655(4) 


16 


0.135720(9) 


6.8 8 


0.135078(8) 


16 


0.135097(5) 


7.4 8 


0.134058(4) 


16 


0.134071(4) 


8.0 8 


0.133168(4) 


16 


0.133173(3) 


9.6 8 


0.131447(3) 


16 


0.131448(2) 


12.0 8 


0.129913(2) 


16 


0.129909(2) 


24.0 8 


0.127261(1) 


16 


0.127258(1) 



In table 1 we list the values of 0, the lattice sizes L/a and the associated 
critical hopping parameters k c at which the numerical simulations have been 
performed. Our final results for the normalization constants are collected in 
table 2. For each value of we quote the number obtained on the larger 
lattice with two errors, the first being the statistical error, which includes the 
uncertainty in the value of k c quoted in table 1. 

The second error is an estimate of the systematic effect which derives 
from the fact that the chosen lattice sizes are not exactly the ones required by 
the normalization conditions. The situation at = 6.0 is exceptional in this 
respect, because L/a = 8 is the correct lattice size and the systematic error 
hence vanishes. The chosen lattice sizes L/a = 12 and L/a = 16 at = 6.2 
and = 6.4 are rather close to the correct ones. In this case the data at 
L/a = 8 may be used to estimate the change in the normalization constants if 
L/a would be lowered to 11 and 14, repectively, which is then quoted as the 
systematic error in table 2. For > 6.8 the error is taken to be the difference 
of the normalization constants calculated at L/a = 16 and L/a = 8. Since 
the effect is of order (a/L) 2 this procedure appears to be safe and presumably 
over-estimates the error. 
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Table 2. Results for Zy and Za 



p 


Zy 


Z A 


6.0 


0.7809(6) 


0.7906(94) 


6.2 


0.7922(4)(9) 


0.8067(76) (23) 


6.4 


0.8032(6)(12) 


0.8273(78)(10) 


6.8 


0.8253(5)(43) 


0.8549(37) (73) 


7.4 


0.8494(3) (34) 


0.8646(20)(48) 


8.0 


0.8667(2) (33) 


0.8812(19)(17) 


9.6 


0.8973(2)(33) 


0.9078(12)(37) 


12.0 


0.9232(2)(26) 


0.9315(11)(16) 


24.0 


0.9656(1)(16) 


0.9692(4) (14) 



Our numerical results are also shown in fig. 4, where we compare them 
with the one-loop expressions [6,23,24], 

Z v = 1 + Z^gl + O(<7 4 ), 4 1} = "0.129430, (6.6) 

Z A = 1 + Z^gl + O(<7 4 ), = -0.116458. (6.7) 

These formulae describe the data rather well for, say, < 0.5, but in the 
range of couplings which is relevant for numerical simulations of physically- 
large lattices this is no longer true. 

The "mean field improved" perturbation expansion [25], 

Z = u {l + (ZW + ±)g 2 P + 0(g*)}, Z = Z y ,Z A , (6.8) 

comes much closer to the data at low values of (3 (see fig. 4). The expansion 
parameter here is Parisi's boosted bare coupling [26], 

9l = 9o/4, (6.9) 

with Uq being the average plaquette in infinite volume at the value of g$ 
considered. We mention in passing that the data are nearly perfectly matched 
by the one-loop formulae (6.6) and (6.7) if we replace g$ by g\. 
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0.9 - 



0.8 - 



0.2 



0.4 



0.6 



0.8 



Fig. 4. Results for Z\ and Za from numerical simulations (filled 
circles) , bare perturbation theory (dotted lines) and "mean field improved" 
perturbation theory (crosses). The solid lines represent the fits (6.10) and 
(6.11). For the numerical data only the statistical errors are displayed. 

In the whole range < go < 1 a good representation of the numerical 
results is given by the rational expressions 



Zv = 



Z A = 



0.7663 gl + 0.0488 g% 
1 - 0.6369 gl ' 

0.8496 gl + 0.0610 g% 



(6.10) 



(6.11) 



1 - 0.7332 gl 

which coincide with the expansions (6.6) and (6.7) to order gl. The fits repro- 
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Fig. 5. Numerical results for the coefficient by. The solid line 
represents the fit (6.12). 

duce the values of Zy and Zx quoted in table 2 with a precision better than 
0.4% and 0.6%, respectively, an exception being the result for Za at [3 = 6.8 
which deviates by 1.35%. Note, however, that the data points shown in fig. 4 
are statistically independent and a statistical fluctuation of this size is hence 
not an unlikely event. 

For future use of our results we suggest to either take the numbers quoted 
in table 2 (where this is possible) or else to employ the fit formulae given above, 
quoting an error of 0.5% for Zy and 1.0% for Z^. These error margins should 
be wide enough to account for all the uncertainties in our calculations. 

6.3 Computation of by 

At non-zero quark mass the renormalized improved currents involve correction 
factors of the form 1 + by^am^ that so far are only known to lowest order 
of perturbation theory [12,17]. The normalization condition (3.13) for the 
vector current, which is also valid for massive quarks, allows us to extract the 
coefficient by by studying the dependence of the ratio fi/ fy(xo) on the quark 
mass m q . In this subsection we report on our results for by obtained on lattices 
of size L/a = 8. At all couplings (3 considered we computed the correlation 
functions at several values of n around n c and extracted by from the linear 
slope of /i//v(T/2) as a function of re -1 . In the range < am q < 0.005 no 
significant curvature was observed in the data. 
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The results of our calculations are plotted in fig. 5. As in the case of the 
coefficients c sw and ca, the data can be represented by a global fit 

l-0.6518g 2 - 0.1226 n ^ 
6v = 1 - 0.8467 g g ' ^ 

which reproduces the tree-level value 6y = 1 a t 50 = 0. 

For studies of pseudoscalar decay constants it would also be desirable to 
know the coefficient 6a- The normalization condition for the axial current 
derived in subsect. 3.3 is however only applicable at zero quark mass. It is 
possible to deduce a more general normalization condition by taking the mass 
term in the PCAC relation into account, but the equation that one obtains 
involves a short distance contribution of order am q which we would not know 
how to separate from the term proportional to 6a- A more sophisticated 
approach is hence required to compute 6a non-perturbatively. We should add, 
however, that a perturbative estimate of 6a may be perfectly satisfactory, if 
one is interested in situations where am q is small (say less than 0.01). 

6.4 Residual cutoff effects in eqs. (3.13) and (3.18) 

Now that the normalization constants Zy and Z^ are known, we can ask 
how large the error term on the right-hand sides of eqs. (3.13) and (3.18) 
is for different choices of the kinematical parameters. In particular, we can 
vary the boundary values of the gauge field, the angles 0k and the insertion 
points xo, yo, and we may also replace one of the quarks created and annihi- 
lated at the boundaries of the lattice by a strange quark with non-zero mass 
(cf. subsect. 3.4). The purpose of such studies is to verify the effectiveness 
of O(a) improvement and also to check that the particular choices made in 
subsect. 6.1 do not lead to uniformly large higher-order cutoff effects in other 
matrix elements of the renormalized currents. 

Several tests, covering all variations mentioned above, have been per- 
formed at (3 = 6.4. For lattice sizes L/a > 8 we found that the cutoff effects 
would amount to changes in the normalization constants no larger than the 
statistical errors quoted in table 2, thus providing another impressive demon- 
stration of the importance and effectiveness of improvement. 
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7. Concluding remarks 



The computation of the current normalization constants and Zy com- 
plements the non-perturbative determination of the 0(a) counterterms in the 
improved action and the improved axial current reported in ref. [18]. In par- 
ticular, physical matrix elements of the axial current in quenched QCD can 
now be obtained with 0(a) improvement fully taken into account and small 
uncertainty in the normalization factor. It should again be emphasized, how- 
ever, that an extrapolation to the continuum limit will always be required, 
even though we have not observed any significant residual cutoff effects in the 
matrix elements considered here. 

The methods we have used in this paper carry over literally to QCD with 
dynamical quarks. As a first step one may be interested in a determination 
of the normalization constants using lattices of size L/a = 8 at all couplings 
§q. The results obtained here suggest that the associated systematic errors 
are quite small (at most 2% in the case of Z&). Simulations of larger lattices 
will however be required for a reliable estimation of the systematic errors and 
for more precise results. 
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DESY-IfH and the staff of the computer centre at Zeuthen for their support. 
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the U.S. Department of Energy (contracts DE-FC05-85ER250000 and DE- 
FG05-92ER40742). Hartmut Wittig acknowledges the support of the Particle 
Physics and Astronomy Research Council through the award of an Advanced 
Fellowship. 
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Appendix A 



We here deduce the Ward identity (3.14) assuming that the axial current 
is conserved (at zero quark mass) and that the operator product expansion is 
valid in a weak sense. 

The axial current conservation amounts to saying that 

(d„A«(x)O) = (A.l) 

for any field product O localized in a region not containing x. It follows from 
this that the integral on the left-hand side of eq. (3.14) is independent of the 
region R (which must contain y and may not have any overlap with the lo- 
calization region of O ext ). We may, for example, take R to be a ball with 
small radius r centred at y. For r — > the integral may then be calculated 
by inserting the operator product expansion of e ahc A a ^{x)A b l/ (y). Up to log- 
arithmic factors the contributions of the composite fields of dimension d are 
proportional to r d ~ 3 . In particular, fields with dimension d > 3 make no con- 
tribution in the limit r — > 0. The only local field with dimension d < 3 and the 
appropriate transformation behaviour under the flavour and space-time sym- 
metries is the vector current V£(y). We thus conclude that eq. (3.14) must be 
valid up to a proportionality constant k. 

To prove that k = 1 we choose R to be the space-time volume between 
two equal-time hyper-planes and integrate over the space components of y. 
For v = the left-hand side of eq. (3.14) is then equal to some matrix element 
of the commutator of the axial charge with itself, while on the other side of 
the equation one has a matrix element of the isospin charge between the same 
states. With the canonical normalization of the charges the matrix elements 
are the same and k must hence be equal to 1. 
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